Vortex glass transition in a random pinning model 
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We study the vortex glass transition in disordered high temperature superconductors using Monte 
Carlo simulations. We use a random pinning model with strong point-correlated quenched disorder, 
a net applied magnetic field, longrange vortex interactions, and periodic boundary conditions. From 
a finite size scaling study of the helicity modulus, the RMS current, and the resistivity, we obtain 
critical exponents at the phase transition. The new exponents differ substantially from those of the 
gauge glass model, but are close to those of the pure three-dimensional XY model. 
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The magnetic field-temperature phase diagram for vor- 
tices in disordered high temperature superconductors has 
been the focus of considerable interest in the past few 
years, centered around e.g. the suggestion of M. P. A. 
Fisher of a possible vortex glass phase with vanishing 
linear resistance [Q. Numerical simulation is often an 
important tool in investigations of phase transitions in 
vortex systems. For example, recent simulations have 
provided valuable information about the phase diagram 
for the case of weak disorder or weak magnetic fields . 
Here evidence for a first order transition separating a 
Bragg glass phase j|, i.e. a dislocation-free solid with 
algebraically decaying translational order, from a vortex 
liquid phase was obtained. However, in the opposite limit 
of strong disorder or strong magnetic fields, the existence 
of a vortex glass phase, i.e. a vortex solid which is topo- 
logically disordered in terms of frozen in dislocations, and 
the critical properties of the vortex glass transition, have 
not been settled. In this paper we study these issues by 
Monte Carlo simulations and finite size scaling. 

Simulations of disordered superconductors have often 
used a so-called gauge glass model B, where the disor- 
der is modeled as a random vector potential added to the 
phase difference of the superconducting order parameter. 
This disorder corresponds to spatially random magnetic 
flux present in the system. Interestingly, recent simula- 
tion results for the gauge glass j| gives similar values 
for the critical exponents as obtained in certain recent 
experiments where a vortex glass has been reported ||. 
However, the gauge glass model has two features that are 
not particularly realistic, namely, it does not assume any 
pinning directly affecting the vortex core energy, and, 
secondly, it is completely isotropic and does not contain 
any net magnetic field. It is at present unknown if these 
details are important for the universality class of the glass 
transition, i.e. if they modify the critical exponents. In 
particular, there is a possibility of anisotropic scaling pro- 
duced by the net field, such that the correlation length 
diverges with different exponents along and perpendicu- 
lar to the field. Some of these and other issues have been 
addressed in the literature. 

An XY model with a random coupling constant and a 



net magnetic field has been simulated using open bound- 
ary conditions ||. However, periodic boundary condi- 
tions are preferable when bulk properties are studied, to 
avoid any influence of the sample surface. Also the effect 
of screening of the vortex interaction has been consid- 
ered. Gauge field fluctuations lead to screening of the 
vortex interaction. This screening is usually rather weak 
in the high temperature superconductors. In models for 
the strong screening limit, simulations give evidence for 
the absence of a vortex glass phase at finite temperatures 
H . A vortex glass scenario without any thermodynamic 
phase transition has also been suggested ||. 

In this paper we consider a random pinning model 
that contains all the pieces necessary to describe the 
static universal critical properties of the vortex glass crit- 
ical point: long-range vortex interactions, strong vor- 
tex pinning, a net applied magnetic field, and periodic 
boundary conditions. The vortex-vortex interaction is a 
full 3D longrange interaction without screening, that be- 
comes applicable when the bare screening length is much 
longer than the vortex spacing, i.e. in the strong field 
limit. The vortex pinning corresponds to uncorrelated 
quenched point disorder, implemented as a position de- 
pendent core energy. This is equivalent to random- T c dis- 
order. The vortex-vortex interaction and the disorder are 
isotropic in the model, but the applied net magnetic field 
breaks the spatial symmetry. We include the possibility 
of anisotropic scaling by allowing for different correlation 
length exponents in the directions parallel and perpen- 
dicular to the magnetic field. The dynamic universality 
class assumed here is that of relaxation dynamics of the 
vortex lines, and we stress that there are more dynamic 
universality classes possible || . 

The random pinning model is denned by the Hamilto- 
nian 

H = \ £ V ( r ~ r ')<l( r ) ■ ^ r ') + \ £ D » ( r M r ) 2 (!) 

r,r' r,/j 

The model is defined on a simple cubic lattice with 
fi = L x L x L z sites, using periodic boundary condi- 
tions in all three directions. The vortex line variables are 
specified by an integer vector field q(r), whose /i = x,y,z 
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component is the vorticity on the link from r to r + e M . 
The partition function is Z = Tr exp(— H/T), where T 
is the temperature, and Tr denotes the sum over all pos- 
sible integers q^, subject to the constraint V • q = on 
all sites, i.e., the vortex lines have no open ends. An 
applied net magnetic field is included as a fixed number 
of vortex lines penetrating the system in the z-direction. 
Three different fillings of the applied magnetic field are 
considered here, i.e., number of vortex lines per link in 
the z-direction: / = 1/2, 1/4, l/ylO. For the irrational 
filling we used the integer number of vortex lines closest 
to fL 2 . The figures below are for / = 1/2. The long 
range vortex-vortex interaction is given by 

nrH "^?E M (2-2cosM (2) 

where K = 4n 2 J (we set J = 1). On each link in 
the system is a short-range point-correlated random pin- 
ning energy with a uniform distribution in the interval 
< D^{y) < K . Hence the lattice constant of the dis- 
cretization lattice in our model corresponds physically to 
a characteristic length scale for variations in the disorder 
energy landscape. 

The Monte Carlo (MC) trial moves are attempts to in- 
sert vortex loops with random orientation on randomly 
selected plaquettes of the lattice. A MC sweep consists 
of one attempt on average to insert a loop on every 
plaquette. The attempts are accepted with probability 
1/(1 + expA_E/T), where AE is the energy change for 
inserting the loop. For half and quarter fillings the ini- 
tial vortex configurations consist of straight lines along 
the z-direction, arranged in a regular lattice in the xy- 
plane. For filling / = 1/vTO the initial configuration 
has straight lines placed at random. For equilibration 
about 10 5 MC sweeps are used, followed by equally many 
sweeps for collecting data. The results were averaged over 
up to 2000 samples of the disorder. Thermal averages are 
denoted by (• • •) and disorder averages by [■■■]. To avoid 
systematic errors in the calculation of squares of expec- 
tation values, two replicas of the system with the same 
disorder are used. 

Superconducting coherence in the vortex line system 
can be detected by calculating the helicity modulus. One 
way of doing this |l(J is to add a term Hq = ^Q 2 to 
the Hamiltonian, where Q M is the total projected area of 
vortex loops added during the simulation. The helicity 
modulus in the direction p is then given by 

y, = i- ^KQD-iQ,) 2 } (3) 

and the RMS current density is defined as J p = 
fj [(Qfj.) 2 ] 1 ^ 2 ■ The linear resistivity, p, is obtained from 
the Kubo formula for the resistance: |[l| 

[w*wo)>] (4) 

t=-t 



where t is MC time, and the voltage is cx AQ^ is 
the net change in the projected vortex loop area during 
a sweep. In the calculation of p the Hq term is not 
included in H. The summation time t is chosen large 
enough that the resistivity is independent of to, but much 
shorter than the length of the simulation. 

We use a generalization of the Fisher-Fisher-Huse ffl 
scaling ansatz to analyze our MC data. At the glass tran- 
sition temperature T c the correlation length in the xy- 
planes, £, and in the z-direction, £ z , and the correlation 
time, t, are assumed to diverge as £ ~ \T— T c \~" , £ 2 ~ 
and t ~ £ z . The anisotropic finite size scaling ansatz |l2] ] 
for the helicity modulus is 

Y x = L*- d -<f 1B (L 1 '''(T-T c ),L t /L<) (5) 
Y z = L 1 - d +<f z (L 1 /»(T-T c ),L z /L<), (6) 

where d = 3 is the spatial dimensionality and / M are 
scaling functions (scaling functions will from now on be 
suppressed in the equations). The current density scales 
as J x ~ L 2 ~ d ~^, J z ~ L l ~ d . The linear resistivity p = 
E/J, where E is the electric field, obeys 

p x ~L d ~*+<-, Pz ~L d - 1 -<-', (7) 

We did a number of tests of equilibration of our simula- 
tions. We followed the standard procedure of calculating 
the "Hamming distance" between two replicas with iden- 
tical disorders for some selected system sizes L, L z . 
We also increased the number of sweeps for equilibration 
with a factor of 10 for a few selected parameter values, 
and obtained no deviations from the data in the figures 
below. The linear resistivity (shown below in Fig. ^) 
gives an estimate of the relaxation time to, which is the 
time where the curves saturate. The values for to gives 
a rough estimate of the required equilibration time. The 
equilibration times used in our simulations are « 8to- 

The first step is to verify that the model has a vor- 
tex glass phase, instead of a Bragg glass phase || that 
is expected for low fillings and weak disorder. Figure 
[j] shows a typical snapshot of a sample configuration for 
T = 0.5(<C T c ), i.e. deep in the glass phase. This calcula- 
tion was done using an exchange MC algorithm |l3] with 
9 uniformly spaced temperatures in [0.5,4.5]. We also 
computed the structure function S(q), and obtain no es- 
sential difference between S in the vortex liquid phase 
and the glass phase, with no indication of a Bragg glass 
phase ||. This demonstrates that the low temperature 
phase of our model, for the filling and disorder strength 
considered, is indeed a glass where the vortex lines are 
frozen in random positions. 

To locate the critical temperature of the vortex glass 
to liquid transition, and determine the critical exponents, 
MC data for the helicity modulus in the x and z direc- 
tions is analyzed by the finite size scaling form in Eqs. 
(§),(§). MC data for different T,L,L Z for / = 1/2 is 
plotted in Fig. @. To determine the anisotropy exponent 
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FIG. 1. Snapshot of a vortex line configuration from our 
MC simulation for / = 1/2 at T = 0.5 (< T c ) which is deep 
in the vortex glass phase. 



£, a sequence of system sizes L z are examined for each 
L. In the relation L z = cL^, both c and £ are varied 
until data curves for and Y z for different system sizes 
L = 6,8, 10, 12 become approximately independent of L 
at a common temperature, which is our estimate for T c . 
Fig. U shows the best crossing obtained in this procedure, 
which gives T c — 4.5 ± 0.1, £ = 1 ± 0.1. The same results 
are obtained from the RMS currents (inset). Equally 
good fits are obtained for c in the interval 1.5 < c < 2. 
The errors are estimated as the interval outside which 
considerably poorer scaling is obtained. 

To determine the correlation length exponent v we 
use fits to Eqs. (H),(H) t° obtain a data collapse for 
system sizes L = 6—12 over an entire tempera- 
ture interval around T c . As a measure of the qual- 
ity of the collapse we define the RMS fit error A = 

El,tJLY^(L,T) - Ux)) 2 ] 1/2 , where x = lM"(T - 

T c ), and / M is estimated by e.g. a cubic polynomial fit 
to the MC data. We did several different types of fits, 
all giving similar results. Figure [| shows the best data 
collapse in a two parameter fit where both T c and v are 
varied independently in the x and z-directions. The best 
fit is obtained for T c = 4.5 ± 0.1, v = 0.7 ± 0.1. A good 
data collapse is obtained, except for the smallest sys- 
tem size, where deviations are obtained below T c . How- 
ever for larger sizes scaling gets better, suggesting that 
the deviation for small system size is a finite size effect. 
The inset shows the RMS fit error as a function of v 



FIG. 2. MC data for the helicity modulus vs. temperature 
for different system sizes L. Inset: RMS current density vs. 
temperature. 

for fixed T c — 4.5. Nearly identical results are obtained 
by instead analyzing data for the RMS currents. The 
same result for v is also obtained from an analysis of the 
derivative of the helicity moduli wrt T. For the lower 
fillings, / = 1/4, 1/ylO, we find similar exponents but 
with larger error bars. 
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FIG. 3. Finite size scaling data collapse of MC data for 
the helicity moduli in the x and z-directions, obtained by 
fitting both T c and v. The data in the ^-direction has been 
shifted to LY X + 1 for clarity. Solid curves are guides to the 
eye. Inset: RMS fit error vs. v at T c — 4.5. 

So far we have presented results for static quantities, 
and we now consider dynamics. The dynamic critical 
exponent is obtained from the linear resistivity in Eq. 
(§, for / = 1/2, T = T c = 4.5, L z = 2L. Fig. | shows 
finite size scaling data collapses according to Eq. (j^) of 
MC data for p x , p z vs. the total MC integration time to 
in the Kubo formula. The inset shows the RMS error 
in a power law fit to the MC data curves, and the best 
fit is obtained for z = 1.5 ± 0.2. This gives a resistivity 
exponent in p ~ t" of s = v{z — 1) w 0.3. 
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FIG. 4. Finite size scaling data collapse of MC data for the 
resistivity, obtained for T c = 4.5, z = 1.45, and L = 8, 10, 12. 
Inset: RMS fit error in the data collapse. 

We will now discuss the values obtained for the critical 
exponents of the random pinning model: £ = 1 ±0.1, v = 
0.7 ± 0.1, z — 1.5 ± 0.2. Our correlation length expo- 
nent v is close to the limiting value at a disordered fixed 
point allowed by the inequality v > 2/d. Within 

the precision obtained, we can not distinguish this expo- 
nent from that of the pure zero field 3D XY model, i.e. 
v « 0.67. Also our dynamic exponent z w 1.5 is consis- 
tent with the zero field 3D XY model with MC dynamics 
for vortex loops jl6]||] . If confirmed, these results suggest 
a scenario where the vortex glass transition in the ran- 
dom pinning model for high fillings and strong disorder 
belongs to the zero field 3D XY universality class. In 
particular, this also suggests that the glass transition is 
driven by similar effectively isotropic, closed vortex loop 
fluctuations as in the zero field case, on top of a glassy 
groundstate. 

Finally we compare our results with other models for 
the vortex glass transition and with experiments. The 
critical exponents obtained here for the random pinning 
model differ considerably from those of the gauge glass 
model g]: v « 1.39, z sa 4.2, s = v(z - 1) » 4.5, and also 
from a random coupling 3D XY model with an applied 
field and open boundary conditions B: v « 2.2, z ~ 
3.3, s w 5.3. Hence, these models appear to belong to 
different universality classes than the random pinning 
model. Experiments on (K,Ba)Bi03 give s w 3.9, and ex- 
periments on untwinncd proton irradiated YB^CuaO^-a 
give s i=a 5.3 H. These experiments show consistency 
with expected vortex glass behavior for tilting the mag- 
netic field away from the c-axis. Unexpectedly, the ex- 
ponents disagree considerably with our results. Further 
work is needed in order to clarify the reasons for this 
discrepancy. 

In summary, we have observed a finite temperature 
continuous vortex glass transition by simulations and 
finite size scaling analysis of a three-dimensional ran- 



dom pinning model. The critical exponents, ( kj 1,1/ w 
0.7, z w 1.5, are surprisingly close to those of the zero- 
field pure 3D XY model, but disagree with the gauge 
glass model and with some experiments. These results 
motivate further theoretical and experimental work. 
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